clc; clear; close all

% 参数设置
b2 = 2*pi;
c  = 11.5;
a_list = linspace(0.33, 5, 500);  % 建议不要从0开始，a=0无动力学意义
n_trans = 500;    % 暂态点数
n_plot  = 100;    % 取用于绘制的点数

x_init = 0.3;
y_init = 0.5;
z_init = 0.6;

% 预分配cell
Xmax_all_cell = cell(length(a_list),1);
a_all_cell    = cell(length(a_list),1);

parfor ai = 1:length(a_list)
    a = a_list(ai);
    x = x_init;
    y = y_init;
    z = z_init;
    
    % 丢弃暂态
    for i = 1:n_trans
        [x, y, z] = SIMM(x, y, z, a, b2, c);
    end

    % 收集分岔点
    X_plot = zeros(1, n_plot);
    for i = 1:n_plot
        [x, y, z] = SIMM(x, y, z, a, b2, c);
        X_plot(i) = x;
    end
    Xmax_all_cell{ai} = X_plot;
    a_all_cell{ai}    = a*ones(1, n_plot);
end

% 合并数据
Xmax_all = cell2mat(Xmax_all_cell');
a_all    = cell2mat(a_all_cell');

% 绘制分岔图
figure
plot(a_all, Xmax_all, 'b.', 'MarkerSize', 2)
xlabel('a')
ylabel('X_{max}')

function [x_new, y_new, z_new] = SIMM(x, y, z, a, b2, c)
    x1 = sin(b2*z) * sin(c/x);
    y1 = sin(b2*x1) * sin(c/y);
    z1 = sin(b2*y1) * sin(c/z);
    x_new = a * x1;
    y_new = a * y1;
    z_new = a * z1;
end
